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Density Induced Phases in Active Nematic 
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We introduce a minimal model for a collection of self-propelled apolar active particles, also called 
as ‘active nematic’, on a two-dimensional substrate and study the order-disorder transition with 
the variation of density. The particles interact with their neighbours within the framework of the 
Lebwohl-Lasher model and move asymmetrically, along their orientation, to unoccupied nearest 
neighbour lattice sites. At a density lower than the equilibrium isotropic-nematic transition density, 
the active nematic shows a first order transition from the isotropic state to a banded state. The 
banded state extends over a range of density, and the scalar order parameter of the system shows a 
plateau like behaviour, similar to that of the magnetic systems. In the large density limit the active 
nematic shows a bistable behaviour between a homogeneous ordered state with global ordering and 
an inhomogeneous mixed state with local ordering. The study of the above phases with density 
variation is scant and gives significant insight of complex behaviours of many biological systems. 

PACS numbers: 87.10.Rt, 05.65.-hb, 64.60.Bd 


Introduction :— Active systems are composed of self- 
propelled particles where each particle extracts energy 
from its surroundings and dissipates it through motion 
towards a direction determined by its orientation. These 
kind of systems are ubiquitous in nature, ranging from 
very small scale systems inside the cell to larger scales m 
[9], vibrated granular media [TOl [TT] etc., and have been 
studied extensively through experiments, theories and 
simulations [HQl]. A collection of head-tail symmetric 
‘apolar’ active particles with an average mutual parallel 
alignment is said to be in a ‘nematic’ state, whereas in an 
‘isotropic’ state particles remain randomly oriented. An 
active system where fluid media do not play important 
role in emergence of ordered state, and thus the hydrody¬ 
namic interactions can be ignored, is called a ‘dry active 
system’ [51 ITOl fTSHIS) . 

Active nature of particles introduces a nonequilibrium 
coupling between density and orientation field, as rep¬ 
resented in terms of curvature coupling current in litera¬ 
ture [HUH HO]. Such coupling in active nematic induces 
unusual properties like large density fluctuation [HEI] 
and growth kinetics faster than 1/3 as in usual conserved 
model [22]. Recent studies of the active nematic found 
a defect-ordered nematic state [231125] as opposed to the 
equilibrium nematic for high particle densities. Recent 
experiment on amolyiod flibrils |26| also found a phase 
with coexisting aligned and disordered fibril domains, 
similar to the defect-ordered state obtained in simula¬ 
tions. But few investigations have been done on the be¬ 
haviours of the active nematic in various density limits, 
especially at low densities. Here we introduce a minimal 
model for two-dimensional active nematic and compare 
various ordering phases of active and equilibrium nematic 
in different density limits. The ordering in the system 
is characterised in terms of a scalar order parameter S 
which is the positive eigen value of nematic order param¬ 
eter Q p in two-dimensions. In the low density limit 
both active and equilibrium systems are in the isotropic 


(I) state with particles randomly oriented throughout the 
whole system (see Fig. Bb) - I), resulting in a small S. 
The Phase diagram of the active nematic as a function 
of packing density C (see Fig. [^a)) shows a jump in 
S close to C = 0.37, whereas in the equilibrium case S 
goes continuously to larger values and an isotropic to ne¬ 
matic (I-N) transition occurs close to C = 0.58. In the 
equilibrium nematic (EN) state particles remain homo¬ 
geneously oriented in the system (see Fig. Bb) - EN). At 
C = 0.37 the active system goes from the isotropic to a 
banded state (BS) where particles cluster and align in the 
perpendicular direction to the long axis of the band (see 
Fig. [^b) - BS-1). With increasing density more number 
of particles participate in band formation (see Eig. Bb) 
- BS-2) and S follows a plateau over a range of density. 
In the large density limit active system shows bistability 
between a homogeneous ordered (HO) (see Eig. Bb)- 
HO) and an inhomogeneous mixed (IM) or local ordered 
state (see Eig. Bb) - IM). This IM state is very similar 
to defect-ordered nematic state in ref. [23H25] . 

Model :— We consider a two dimensional square lat¬ 
tice. At each vertex ‘T we define an occupation variable 
Ui, which can take values I (occupied) or 0 (unoccupied), 
and an orientation variable Oi , which lies between 0 and tt 
because of the apolar nature of the particles. Each parti¬ 
cle interacts with its nearest neighbours through modified 
Lebwohl - Lasher Hamiltonian [28] 

V, = —e ^ UiUj cos2{0i — Oj) (I) 

<ij> 

where e is the interaction strength between two neigh¬ 
bouring particles. This model is analogous to the diluted 
XY-model with nonmagnetic impurities [29], where im¬ 
purities and spins are analogous to vacancies and parti¬ 
cles respectively in the present model. 

Orientation evolves through Monte - Carlo (MC) up¬ 
dates [30] following the Hamiltonian in Eq. Unlike 
the diluted XY-model, particles also move on the lat- 
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FIG. 1: (Color online) (a) Plot of scalar order parameter S vs. packing density C for active (circles and triangles) and equilibrium 
(continuous line) nematic for system size 512 x 512. Equilibrium system goes smoothly from isotropic (I) to nematic (EN) 
state. Active system goes from isotropic (I) to banded state (BS) (small jump in S) followed by either an inhomogeneous 
mixed (IM) (triangles) or a homogeneous ordered (HO) (circles) state, (b) Snapshots of particle inclination towards the 
horizontal direction. Color bar ranging from zero to one indicates parallel and perpendicular inclinations respectively towards 
the horizontal direction. White regions signify unoccupied sites. (I) is isotropic state at low density (C = 0.36), (BS-1) 
(C = 0.38) and (BS-2) (C = 0.56) are two banded state configurations, (IM) is inhomogeneous mixed, (HO) is homogeneous 
ordered and (EN) is equilibrium nematic state at high density (C = 0.76). 


tice. Depending on the type of movement we define two 
kinds of models, (i) ‘Equilibrium model’ (EM) - a par¬ 
ticle can diffuse to any unoccupied nearest-neighbouring 
site, and therefore satisfies the detailed balance condi¬ 
tion. (ii) ‘Active model’ (AM) - a particle can move to 
only those unoccupied nearest-neighbouring sites which 
are in the direction that makes the least inclination with 
the particle orientation. Details of the model and particle 
movement are shown in Supplemental Material m sec¬ 
tion 1. The asymmetric move of the active particles does 
not staisfy the detailed balance condition and arises in 
general because of the self-propelled nature of the parti¬ 
cles in many biological [151 E2] and granular systems m- 
These moves produce an active curvature coupling cur¬ 
rent in coarse-grained hydrodynamic equations of motion 

naEo]. 

Numerical details :— We consider a collection of N 
particles with random orientation Oi G [0,7r], homoge¬ 
neously distributed on a L x L square lattice {L = 
150, 256, 512) with periodic boundary condition. The 
packing density of the system is C = N/{L x L). We 
choose a particle randomly and move it to an unoccupied 
neighbouring site, followed by an orientation updation 
through MC. We use 10^ MC steps to evolve the system 
to its steady state and all the results have been obtained 
by averaging over next 2 x 10^ MC steps. Twenty four 
realizations have been used for better statistics. 

We calculate the scalar order parameter 


5 ' = 



rii cos(26>i))2 + sin(26>i))2 (2) 

i 


which is small in the isotropic state and close to 1 in the 
ordered state. First we calculate S for EM as a function 
of inverse temperature (3 = l/ksT for different densities. 
As shown in Supplemental Material m section II, the 
critical temperature Tc is approximated as Tc{S = 0.4). 
Critical temperature Tc{C) decreases with the lowering 
of the packing density C , similar trend is found in the 
study of diluted XY-model [29] for varying nonmagnetic 
site density. In rest of our calculations temperature is 
kept fixed at /3e = 2.0 and packing density C is varied 
from small values to complete filling C = 1.0. 

Phase diagram :— At low densities, C < 0.37, the 
active system is in the isotropic state where the parti¬ 
cles with random orientation remain homogeneously dis¬ 
tributed throughout the system, and therefore S holds 
vanishingly small values. The jump occurs in S at 
C = 0.37. For C > 0.37 particles cluster in, and both or¬ 
dered state with high local density and disordered state 
with low local density coexist (see Fig. Bb) - BS-1). 
Mean alignment inside the band is perpendicular to the 
long axis of the band. In the moderate density lim¬ 
its, band formation in more favourable than lane for¬ 
mation (mean alignment parallel to the long axis of the 
structure) because the number of particles that can have 
translational motion is much larger in the banded state, 
and therefore entropy favours band formation. Similar 
mechanism create the bending instability close to order- 
disorder transition in the active polar systems [33l|34|. 
The above transition from I state to BS occurs at density 
lower than the corresponding equilibrium I-N transition 
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density C ^ 0.58 (see Fig. [^a)). 

As we further increase C, unlike the equilibrium sys¬ 
tem where S increases monotonically with C, the active 
system shows very small change in S for a range of den¬ 
sity. This plateau like appearance of S with variation 
in C is very similar to plateau phase in magnetization 
versus field curve of magnetic systems [35]. If an energy 
gap exists between two consecutive magnetic states, a fi¬ 
nite field is required for the magnetic system to go from 
the lower to the higher state. So until that finite field 
is applied, the increasing field keeps the system magne¬ 
tization to be unchanged, and the system is called to be 
in the plateau phase. With increasing packing density in 
the plateau regime of the active nematic more particles 
participate in band formation (see Fig. Bb) - BS -1 and 
BS- 2 ). On further increment of density, close to equi¬ 
librium I-N transition C ^ 0.58, transverse fluctuations 
lead the system to a mixed state mm- 

In the large C limit active system shows a bistable 
behaviour with two distinct steady states; first, a state 
where S is large and real space configuration is ‘homoge¬ 
neous ordered’ (HO), and the second, an ‘inhomogeneous 
mixed’ (IM) state where S realizes some moderate val¬ 
ues. In the HO state though the particle orientation is 
homogeneous, large density inhomogeneity exists in the 
system (see Fig. Bb) - HO). IM state is a local ordered 
state with many aligned clusters of high particle den¬ 
sity. The system consists of many such aligned clusters 
of high density separated from low density disordered re¬ 
gions and mean alignment in each cluster is in different 
directions (see Fig. Bb) - IM). IM state is similar to 
the defect-ordered state recently found in the study of 
ref. |23ll25] . with large number of ±1/2 defects in high 
density active nematic. 

Renormalised mean field study for small S: — We also 
calculate the jump in the scalar order parameter S and 
the shift in the transition density using the Renormalised 
mean field (RMF) method of the coupled coarse-grained 
hydrodynamic equations of motion for the number den¬ 
sity p(r,t) = “ Hc^(t)) and the order parame- 

ter Wij{r,t) = p{r,t)Q{r,t) = - ^Sij)S{r - 

for active nematic [HUSO]. 


dtp = aoVjVjWjj + DpV'^p (3) 


and 



FIG. 2: (Color online) g 2 {r) vs. r on log-log scale at different 
densities, (a) Active nematic: (Q) and (□) at low densities 
g 2 {r) decays exponentially, (o) and (±) at intermediate den¬ 
sity g 2 {r) decays algebraically, (A) homogeneous ordered (al¬ 
gebraic decay), and inhomogeneous mixed (abrupt decay) at 
high density, (b) Equilibrium nematic (Q) and (□) exponen¬ 
tial decay of g 2 {r) at low densities and (o) and (±) algebraic 
decay of g 2 (r) at high densities. 


m or from microscopic rule based model [20|. Den¬ 
sity equation is a continuity equation dtp = —V • J, 
because the total number of particles is conserved. Cur¬ 
rent Ji = — aoVjw^j — DpVip^ where the first term con¬ 
sists of two parts, an active curvature coupling current 
Ja oc aopS/jQij and anisotropic diffusion Jpi oc QijVip^ 
which can also be present in the equilibrium model. The 
second term in density equation is an isotropic diffusion 
Jp2 oc Vp term. First two terms in the order param¬ 
eter equation is the alignment term. We choose 
<^i(p) = — 1 ) as a function of density which changes 

sign at some critical density pijsf. Third term is coupling 
to density and last term is diffusion in order parameter 
and written for equal elastic constant approximation for 
two-dimensional nematic. 

A homogeneous steady state solution of Eqs. and 
gives a mean field transition from isotropic to nematic 
state at density piN where ai{p) changes sign. Us¬ 
ing RMF method we calculate the effective free energy 
feff{S) close to order-disorder transition when S is small. 
We consider density fluctuations 5p and neglect order pa¬ 
rameter fluctuations. The effective free energy 


dt^ij = {<^1 (p) - 0^2 (w : ^)}^ij 

+ P ~ ^ ^ (4) 

where, nic, = {cos{0a),sm{0a)) is the unit vector along 
the orientation and Ra(t) is the position of parti¬ 
cle a. We can obtain the number density p by coarse- 
graining C over small sub volume. Eqs. and can 
be obtained either from symmetry arguments as in ref. 


= +jS* (5) 

where 62 = <ai(p) ± a'i{p)c^ where c is a constant. 
(a'i(p) = dai/dp\p^, 63 = ^4 = ^<^ 2 - Both 

63 and 64 are positive. A detail calculation for /e// is 
shown in Supplemental Material [31] section HI. The den¬ 
sity flcutuations introduce a new cubic order term pro- 
tortional to the activity strength ao, in the free energy 
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feff{S). The presence of such term produces a jump 

A5' = 5'c = at a lower density pc = PinO- ~ i^)* 
This type of jump and shift in transition because of fl- 
cutuations are also called as fluctuation dominated first 
order phase transition in statistical mechanics m and 
widely studied in many systems [38]. The jump in S and 
the shift in pc is proportional to the activity parameter 
ao and for uq = 0 we recover the equilibrium transition. 

Two-point orientation correlation function 
:— To further characterise the system we also 
calculate the two-point orientation correlation 
g2{r) =< UiUi^r cos [2 {Oi - Oi^r)]/^inini > at 

different packing densities, where < . > signifies an 
average over many realisations. In Fig. [^we plot ^ 2 (^) 
vs. r on log-log scale, for C = 0.30, 0.36, 0.38, 0.52 and 
0.82 for active model and C = 0.30, 0.56, 0.58 and 0.78 
for equilibrium model. For very small packing density 
C < 0.37, g 2 {r) decays exponentially in the active 
systems. Therefore the system is in short-range-ordered 
(SRO) isotropic state. At C = 0.38, ^ 2 (^) decays fol¬ 
lowing the power law ^ 2 (^) — and therefore the 

system is in a quasi-long-range-ordered (QLRO) state. 
At high packing densities correlation functions confirm 
the bistability in the active systems. For C = 0.82 
(see Fig. [^a)) ^ 2 (^) shows power law decay in HO 
state, whereas in IM state ^ 2 (^) decays abruptly after 
a distance r. The abrupt change in g2{r) at a certain 
distance indicates the presence of local ordered clusters. 
In contrast, the equilibrium systems show a transition 
from SRO isotropic state at low density C ^ 0.56 to 
QLRO nematic state at C ^ 0.58, similar to Berezinskii 
- Kosterlitz - Thouless (BKT) transition |39l [40] in the 
diluted XY-model m- 

Orientation distribution :— We also compare the 
steady state properties of active and equilibrium mod¬ 
els in the high density limit. First we calculate the 
steady state (static) orientation distribution P{0) from 
one snapshot of particle orientation. Both active HO and 
equilibrium nematic show a Gaussian distribution of ori¬ 
entation (see Fig. [^a)). Hence in HO state orientation 
fluctuations of particles are of same kind as in equilib¬ 
rium model and the system is in QLRO state. Distribu¬ 
tion P{0) in the IM state is very broad and spanning the 
whole range of orientation. Hence the system has many 
local ordered regions with different orientations. 

We also calculate the time averaged distribution of 
mean orientation of all the particles P{0m) in active HO 
and equilibrium nematic states. PiOm) is calculated from 
mean of all particle orientaions averaging over a long time 
(from 10^ to 3 X 10^) in the steady state. P{0m) for HO is 
narrow in comparison to the broad distribution for equi¬ 
librium model (see Fig. ib)). Narrow distribution of 
P{0m) implies that orientation autocorrelation in active 
system decay slowly as comapared to the corresponding 
equilibrium model which is in agreement with the slow 
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FIG. 3: (Color online) (a) Steady state orientation distribu¬ 
tion P{0) of particles for HO, IM (active) and EN (equilib¬ 
rium) states at high density. Lines are fit to Gaussian dis¬ 
tribution for both HO and EN states. IM state shows very 
broad distribution of P{0). (b) Plot of mean orientation dis¬ 
tribution P{0m) averaged over a long time in the steady state 
for HO (active) and EN (equilibrium) states. P{0m) is very 
broad for EN in comparison to HO. 


decay of velocity autocorrelation in bacterial suspension 

SI]- 

Summary :— In this letter we have introduced a min¬ 
imal model for active nematic and found three distinct 
phases with the variation in density. At low densities 
the active nematic is in disordered isotropic state with 
very small correlation between the particles. With in¬ 
creasing density active nematic undergoes a fluctuation 
induced first order phase transition from the isotropic to 
the banded state where large number of particles partic¬ 
ipate in band formation. Large density fluctuations in 
the active systems change the nature of the transition 
and shift the transition density to smaller value as com¬ 
pared to the equilibrium isotropic nematic transition. At 
large densities equilibrium nematic is in QLRO nematic 
state, whereas active nematic goes from the banded state 
to either the homogeneous ordered (high S) or the inho¬ 
mogeneous mixed (moderate S) state. This inhomoge¬ 
neous mixed state is similar to the phase with coexist¬ 
ing aligned and disordered fibril domains found in recent 
experiment m- Experiments on thin layer of agitated 
granular rods, elongated living cells, bacterial colony of 
apolar B. suhtilis etc. at different densities can realize 
the different phases we found here. In the present model 
we have frozen the motion of the active particles in the 
transverse direction, i.e. the activity strength is kept 
large. It will be interesting to see the evolution of differ¬ 
ent phases with the particles having a small probability 
to move in transverse directions as well. 
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Supplemental Material 


I. MODEL FIGURE 
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FIG. 4: (a) Two dimensional square lattice with occupied (n = 1) or unoccupied (n = 0) sites. Filled circles signify the 
occupancy of respective sites. Inclination of the rods towards the horizontal direction show respective particle orientation 
Oi G [0,7r]. (b) Equilibrium move : particle can move to any of four neighbouring sites with equal probability 1/4, (c) Active 
move: particle can move to either of its two neighbouring sites with probability 1/2, if unoccupied, in the direction it is more 
inclined to. 
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II. ESTIMATE OF CRITICAL TEMPRATURE Tc(C) IN EQUILIBRIUM MODEL 



Pe 

FIG. 5: Plot of S vs. inverse temperature (Se for different densities C for equilibrium model. System goes from isotropic (small 
S) to nematic (large S) state. Vertical dotted line shows the variation in S for fixed /3e = 2.0 at different densities. Crtical 
temperature is approximated as Tc{S = 0.4). Inset: change in Tc as a function of density C. 


III. RENORMALISED MEAN FIELD (RMF) STUDY OF ACTIVE NEMATIC FOR SMALL SCALAR 

ORDER PARAMETER S 


In this section we will write an effective renormalised mean field free energy for scalar order parameter S for small S. 
We keep the density fluctuations and ignore the order parameter fluctuations in the coupled hydrodynamic equations 
of motion for active nematic. Density fluctuation produces a cubic order term in the effective free energy for scalar 
order parameter S and such term produces a jump in S' at a new transition density pc lower than equilibrium I-N 
transition point p/w- Shift in transition density and jump AS is directly proportional to the activity parameter uq 
and we recover equilibrium limit for zero ao- 

We write coupled hydrodynamic equations of motion for density p and order parameter w = pQ where nematic order 
parameter [T] 


O^r f\=^( cos26>(r,i) sm2d{r,t) \ 
^ ^ ysin2^(r, t) — cos2^(r,t) J 

0 being the coarse grained angle at position r and time t. Density equation 

dtp = aoViVjWij + DpV'^p 


(6) 

(7) 


and order parameter equation w 

= {«! {p) - Q ;2 (w : w)} Wij + /3 ^ViVj - P + -DwV^w^ (8) 

Density Eq. [^is a continuity equation dpjdt = —V • J, where J has two parts, active and diffusive. Details of these 
two currents are given in the main text, ao is the activity parameter, present because of self-propelled nature of 
the particles, p is the coupling of density in w equation. Dp and are the diffusion coefficients in density and 
order parameter equations respectively, ai{p) and 0^2 are the alignment terms and ingeneral depends on the model 
parameters. For metric distance interacting models [2] ai{p) is a function of density and changes sign at some critical 
density. We choose ai{p) = — 1 and 0^2 = 1. Steady state solution of homogeneous Eq. j^is p = po? we add small 

perturbation to mean density p = p^ ^ 5p. In the staedy state density fluctuation 5p can be obtained from Eq. 

aoViVjWij + DpV‘^6p = 0 

^ aoS/jWij + DpVidp = constant = Ci 


( 9 ) 









( 10 ) 


where wn = —W 22 = f cos(26>) and W 12 = W 21 = f sm{20) and keep the lowest order terms in S and 0 

djp = -^d,S ^ 5p{x) = -^S + c 


and 


dy5p=^dyS^5p{y) = ^S + c^ 


ao 


( 11 ) 


Here we assume nematic is aligned along one direction and there is variation only along the perpendicular direction. 
Hence we can choose either of equations or Two constants c and ci are the value of density where nematic 
order parameter is zero. 

We use Eq. 10 and substitute the solution for density in equation for and obtain an effective equation for S. 


dtS = { ai (p) - ^a2S'^ \ S + 0{V^S) + 0{V^p) 


( 12 ) 


We neglect all the derivative terms and keep only polynomial in S, i.e. we neglect higher order fluctuations. We 
can do taylor expansion of ai{p) about mean density po, = <^i(po + ^p) = <^i(Po) + Ci[6p^ where a[ = ^^|po- 

Substitute the expression for 6p from Eq. hence 


dtS = 


|ai (po) + a'l^p - 




We can write an effective free energy for S 


hence 


dtS=- 


SfeffjS) 

6S 




Therefore 


/e//(5) = -|52 






(13) 


(14) 


(15) 


(16) 


where 62 = (po) + 63 = and 64 = ^ 0^2 and c is a conatant. Hence fluctuation in density introduces a 

cubic order term in effective free energy /e//(5'). Effective free energy in Eq. 16 is similar to Landau free energy 


with cubic order term [3]. We calculate jump AS' and new critical density from coexistence condition for free energy. 
Steady state solutions of order parameter (S = 0 for isotropic and S 7 ^ 0 for ordered state) are given by 


^-0 = (-62 - bsS + 6452) 5 = 0 

Non-zero S is given by —62 — b^Sc + = 0 . Coexistence condition implies 

fefPSc) = 0 - jS, + ^52) 52 = f,fPS = 0) = 0 

hence we get the solution 


Sc = 


% 

bs 


(17) 


(18) 


(19) 


and 


b 


C _ 

2 — 


24 

964 


( 20 ) 
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Hence the jump at new critical point is AS' = Since 64 > 0 and hence 62 < 0, the new critical density 


2 ^ 2 ' 

Pc = Pin ( ^ ^ ^ 


( 21 ) 


is shifted to lower density in comparison to equilibrium pijsf. Eq. gives the expression for new transition density 
as given in main text. Hence using renormalised mean field theory we find a jump AS at a lower density as compared 
to equilibrium I-N transition density. 
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